Noname manuscript No. 

(will be inserted by the editor) 



o 

(N 



A slip model for micro/nano gas flows induced by body 
forces 

Q.D. To • C. Bercegeay • G. Lauriat • C. Leonard • G. Bonnet • 



Received: date / Accepted: date 



(N 



V 

i 

q 

O 

^ ■ 

Oh. 



^, 

O 
O 



X 



Abstract A slip model for gas flows in micro/nano 
channels induced by external body forces is derived 
based on Maxwell's collision theory between gas molecules 
and the wall. The model modifies the relationship be- 
tween slip velocity and velocity gradient at the walls 
by introducing a new parameter in addition to the clas- 
sic Tangential Momentum Accommodation Coefficient. 
Three-dimensional Molecular Dynamics simulations of 
helium gas flows under uniform body force field between 
copper flat walls with different channel height are used 
to valid the model and to determine this new parame- 
ter. 

Keywords Rarefied effect • Kinetic Maxwell model • 
External volume force • Slip model • Tangential Mo- 
mentum Accommodation CoefHcient • MD simulation 



1 Introduction 

The velocity of a fluid close to a solid wall is always 
different from the wall velocity even if the latter is 
perfectly diffusive. Especially, when the channel height 
is decreased to that of MEMS or NEMS devices (Mi- 
cro/Nano Electro- Mechanical Systems) , this phenomenon 
becomes highly important and must be taken into ac- 
count. The Knudsen number Kn, the ratio between the 
mean free path A and the characteristic length of the 
channel H, is the relevant parameter to quantify the 
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NOMENCLATURE 



Mean free path 

Channel height, width and length 

Knudsen number 

Number density, mass density 

Molecular diameter 

Slip length, dimensionless slip length 

Cartesian coordinate 

Normalized coordinate 

Tangential Momentum Accommodation Coefficient 

Tangential velocity, normalized tangential velocity 

Slip velocity, normalized slip velocity 

Reference velocity 

Molecules going upward A'^+ and downward N~ 

with respect to the control surface s 

Total number of molecules passing through s 

Average velocity of molecules going upward < dJ > 

and downward < v^ > with respect to s 

Gas velocity near the wall, velocity of the wall 

Mean collision time, thermal speed 

Boltzmann constant, absolute temperature 

Molecular mass, acceleration along x-axis 

Slip parameters of the present model 

Gas velocity at distance A from the wall 

Gas viscosity, scaled viscosity, kinetic theoretical 

viscosity 

Potential energy of atom i 

Potential functions of Embedded Atom Model 

Parameters of Lennard Jones potential between 

molecules a and b 

Distance between two molecules i and j. 



slip effects. The mean free path A is usually defined 
as the average distance that molecules travel between 
collisions and equal to 
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where n is number density and d is the effective molecu- 
lar diameter. According to Maxwell's model [T], the slip 



length Ls in continuous fluid mechanics can be deter- 
mined via the Tangential Momentum Accommodation 
Coefficient, also denoted by TMAC or a^ as follows 

L, = ^^^A (2) 



is calculated by the formula 



Vsli 



(3) 



The slip velocity, VsUp 

dv 
oz 

The term ^ is the normal derivatives of the tangen- 
tial velocity component at the walls, assuming that the 
normal to the wall is in the z-direction. If the velocity 
and coordinate in Eq. ^ are scaled with a reference 
velocity Vref and the channel height H, we have 
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where VsUp = VsUp/vref and z = zjH . The term L 
called the dimensionless slip length and equal to 

LI = ^^^Kn 



IS 



(5) 



when the Maxwell model is used. Consequently, for a 
given accommodation coefficient, Eq. ([5]) predicts that 
L* is proportional to Kn. When Kn tends towards zero, 
we recover the no slip condition and when Kn increases, 
the slip effect increases. The Maxwell model is widely 
used to describe the slip at the walls because it only 
needs one parameter only, (t„. 

The physical meaning of ct„ in (J2l5p is that if M molecules 
arrive at the wall, a^M of them are reflected diffusively 
and the remaining (1 — a^^M molecules are reflected 
specularly. Based on Eq. ([2]), one can determine TMAC 
by either experiments or Molecular Dynamics (MD) 
simulations. Some experimentalists controlled macro- 
scopic quantities such as pressure and mass flow rate 
and made use of the relationship with the slip velocity 
to find TMAC (see e.g. 013]). Arkihc et al. studied 
flows of nitrogen, argon and carbon dioxide through 
rectangular silicon channels and found TMAC ranging 
between 0.75 and 0.85. Colins et al. [3] worked on the 
couples silicon-nitrogen and silicon-helium and found 
a relatively high coefficient, 0.93. Recently, Maali and 
Bhushan [3] studied confined air flow between a spher- 
ical glass particle glued to the cantilever of an atomic 
force microscopy (AFM) and a glass plate. They let 
the cantilever oscillating and measured the hydrody- 
namic damping factor. Consequently, a value of the slip 
length of about 118 nm was obtained. It corresponds 
to an accommodation coefficient of 0.72. Direct mea- 
surements of TMAC on the couple He-Cu by Seidl (see 
[51[5]) stated that the coefficient depends on the colli- 
sion angle between gas and solid, ranging from 0.6 to 



1.0. On the other hand, Cao et al. ^^^ using MD 
approaches to simulate flows, revealed that TMAC at 
Ar-Pt interface can be as small as 0.2 and influenced 
by temperature and surface roughness. Arya et al. |10| 
simulated directly the collisions between gas and wall. 
Their results showed the dependence of TMAC on the 
wall's lattice structure. It remains constant as long as 
the drift velocity is small enough (less than 100 m/s). 
Finger et al. [5] used similar method as [TO] to show 
that TMAC is affected by the adsorbed layer and their 
results matched the previous experiment of Seidl on the 
couple copper-helium. Generally speaking, both exper- 
imental and MD works gave a rather scattering results 
of TMAC and reflected the dependency of TMAC on 
many factors. Hence, this coefficient should be under- 
stood as an effective one and to be used with caution. 

Extensions of Maxwell's model to different configura- 
tions have been considered in the past. When the Knud- 
sen number increases, Beskok (see [lllf^ and the ref- 
erences cited therein) argued that higher order of the 
Knudsen number and higher order derivatives of veloc- 
ity must be used in the slip equations. For curved sur- 
face, Lockerby et al. [T^] accounted for the contribution 
of normal velocity in the viscous stress. In all the afore- 
mentioned works, the gas molecules are not subject to 
any external forces before arriving at the wall, which is 
not applicable in the presence of volume force field such 
as gravity, electrostatics, etc... With these force fields, 
the slip velocity is different from fiows driven by pres- 
sure gradient and cannot be simply described by ^ . In 
what follows, on the basis of kinetic theory, we derive 
the new slip equation that accounts for the external 
body force field. 
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Fig. 1: Collision model between gas molecules and a 
solid wall. N^ and N^ are respectively the number 
of molecules going downward and upward and passing 
through the surface s near the wall within one unit 
time. If there is no accumulation of gas near the wall, 
N- = iV+. 



2 Slip model for micro/nano gas flows induced 
by body force 

In the following derivation, the gas flow is assumed to be 
isothermal so that the influence of temperature on the 
slip velocity is not taken into account. Using the similar 
approach of Beskok (see [H]), let us consider a control 
surface s near and parallel to the wall (see Fig. [1]) . For 
a unit time, there are N gas molecules passing through 
the surface which are composed of N~ molecules going 
downward and iV+ going upward with respectively av- 
erage tangential velocities, < v~ > and < v^ >. The 
gas average tangential velocity at the wall Vs may be 
written as 



following reasons 

- the wall in the model is idealized as a surface and 
the gas wall collisions only take place at this surface. In 
fact, the wall also has an atomistic structure and the 
interaction force must be taken into account at distance 
of several molecular diameters. 

- after arriving at the wall, the gas molecules can stay 
near the wall during a flnite duration of time before 
leaving. 

The value of /3 is expected to be close to unity. The 
slip equation with the corrected term becomes 
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N+ <vt > +N~ < v~ >= Nvs 
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The molecules that go upward are those that previously 
went downward and were reflected at the surface s. Be- 
cause the reflection is either diffusive or specular by the 
fractions cr-„ and 1 — cr^ , the following relation holds for 
a wall moving at velocity v.^: 

N+ < v+ >= (1 - <7^)N- <v^ > +a.^N~v^ (7) 

which is equivalent to 

Nvs = N~ [< v^ > +(1 - ay) < v^ > +ayVy,] (8) 

Without external volume force, the A^^ molecules col- 
liding with the wall come from one mean free path A 
away from the wall without collision so that their ve- 
locity does not change < v~ >= v\. We also assume 
that N- =N+ = N/2. It follows 

2vs = [(2 - ay)vx + ayVyj] (9) 

The Taylor development of v\ near the wall, vx = v^ + 

dv 



Equation P^ is the new slip equation for general flows 
in which the body force is involved. In what follows, we 
study a particular case where the gas flow of viscosity /i 
is induced by a constant body force pj^ along one direc- 
tion X. Without pressure gradient, the velocity profile 
is given by 
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which after combined with the new slip model equation 
(fT2|) yields the dimensionless form 
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In the above equation, jj* is the scaled viscosity and 
fl is the kinetic theoretical viscosity [THlfH] defined as 
Jl = \p\c. The dimensionless slip length L* can be 
calculated accordingly 



A 1^1^, gives the relation for slip velocity VsUp defined L* = aKn [/3yU*Kn + 1] with L* = Ls/H 



as Vs 



(first order Maxwell's relationship) 



(15) 
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In the presence of uniform volume force, i.e. in the case 
where a constant acceleration ^x is applied on each gas 
molecule, < v~ > is no longer equal to v\. When im- 
pinging at the wall surface, the term ^^P < t > should 
be added to the average tangential velocity 

<v- >=vx+-fx/3 <T> . (11) 

Above, < T > is the mean time for the molecules to 
arrive at the surface after the previous collision and is 
assumed to be A/c where c is the thermal speed of the 
gas. For gases in local equilibrium, the thermal speed 
is estimated by c = y^2kBT/m. The constant /?, intro- 
duced in the last term of (jll|) . can be seen as a factor 
which accounts for the differences between the idealized 
conditions used to derive the slip model and the realis- 
tic ones. In reality, these differences can be due to the 



We can deduce that the derived slip length is second- 
order dependent of the Knudsen number. The influence 
of the volume force on the slip length becomes thus im- 
portant when Kn is large enough. The ratio L*/Kn is 
no longer equal to the constant a but is dependent on 
the channel characteristic dimension H via a composite 
parameter (/i*Kn)^^. If Kn is increased and the vari- 
ation of fluid viscosity fj, is small, we shall observe an 
increase in the ratio L*/Kn. In what follows, we use the 
molecular dynamics approach to valid this prediction 
and determine as well the two new model parameters 
a and (3. 



3 Molecular dynamics validation 

In order to recover the dependence of slip velocities 
on the Knudsen number and volume force, we simulate 
a flow induced by uniform external volume force. The 
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Fig. 2: One dimensional flow induced by external vol- 
ume force. The simulation domain is periodic in direc- 
tions x^y. The external volume force in direction x is 
produced by applying additional constant acceleration 
7a; on each gas atoms. 



gas/solid couple investigated is helium (He) and copper 
(Cu). The choice of the couple He-Cu is motivated by 
two reasons: firstly, the accommodation coefficients de- 
rived experimentally and numerically are shown to be in 
good agreement (see for example ref. [S]) and, secondly 
the interaction potentials between these atomic species 
have been widely used (see ref. [TB,16,17,5T8i). How- 
ever, the present method can be applied to any gas/wall 
system as long as the interaction potential between the 
atoms is properly provided. Because the slip length de- 
pends on what happens at the interface, the final results 
are strongly affected by the choice of the system. 

The channel is made of two parallel solid walls, each of 
them contains 9 layers (or 4 lattice units) of Cu atoms 
placed at fee lattice sites. The lattice constant is cho- 
sen initially equal to 3.615 A. The interaction forces 
between Cu atoms are derived from the EAM potential 

USES] 
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where Vi is the potential energy of atom i composed 
of the binary potential 4>{rij) and embedded potential 
F accounting for electron density contribution p^. The 
interaction forces for the couples He-He and He-Cu are 
derived from the widely used Lennard- Jones potential 
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The parameters used in this paper are those from [THl 

0, EHe-He = 0.00088 cV, aHe-He = 2.28 A, EHe-Cu = 

0.0225 eV, ane-Cu = 2.29 A. The He-Cu parameters 



Table 1: Size of the simulation box and number of he- 
lium and copper atoms for different Knudsen numbers 

Kn H [A] L [A] B [A] Cu atoms He atoms 



0.046 361 
0.064 260 
0.098 174 



361 
260 

174 



180 
130 

87 



230 700 
51 840 
23 040 



58 800 
22 680 
6 624 



have been calculated using the Lorentz-Berthelot mix- 
ing laws from the He- He and Cu-Cu parameters in |18j . 
The last atom layer at the two walls is fixed (see Fig. 
[2). The remaining atoms of the wall and the gas are 
maintained at the same temperature 120 K by a usual 
scaling method after removal of the mean velocity. On 
average, the wall is kept immobile, i.e v^j = 0. 

In this work, the global number density of the gas are 
kept constant n = 0.0026 A'^ (or p == 1.78 x 10" ^^^ 
pg/A"^) while the channel height H and acceleration 
7a; is varied to obtain results for different global Knud- 
sen numbers and volume forces /. At the small den- 
sity number and temperature of our simulations (e.g. 
120 K), the helium is in gaseous phase. Three values 
oi H ~ 361, 260 and 174 A corresponding respec- 
tively to Kn=0.046, 0.064, 0.098 are considered. These 
three Kn numbers are chosen to fall into the range 
(0.01,0.1) in order to assure the slip flow regimes in 
our test cases according to [TT]. The two other dimen- 
sions of the simulation box along x-axis (the flow direc- 
tion) and along y-axis denoted shortly by length L and 
width B are L = H and B = H/2. For each channel 
height, the acceleration applied on each atom is var- 
ied as 7a; = 0.0036, 0.012, 0.024 A/ps^. The computa- 
tions are carried out by using LAMMPS, an open source 
parallelized code [T^ on an IBM Power6 machine. The 
equations of the particle motion are integrated using 
a Leapfrog- Verlet algorithm with a time step 0.002 ps. 
The steady state is achieved after 5 x 10^ time steps 
and it takes another 2 x 10® time steps for the average 
process. All the models are constructed in 3D with the 
parameters given in Table [1] Due to the high density of 
the solid wall, an important number of Cu atoms are 
considered. The largest model, case Kn=0.046, involves 
289 500 molecules and takes 8 hours of computation on 
512 processors. 

In order to obtain accurate solutions, the channel has 
been divided in a sufficiently large number of layers 
for the average procedure. For 7a; = 0.012 A/ps^ and 
different Knudsen numbers, the velocity profiles in the 
upper half channel are plotted in Fig. |31 We find good 
agreement between analytical solution (Eq. [T5| and the 
numerical results in the major part of the channel. The 




Kn = 0.046, Numerical results 
-Kn = 0.046, Curve fitting 

Kn = 0.064, Numerical results 
-Kn = 0.064, Curve fitting 

Kn = 0.098, Numerical results 
-Kn = 0.098, Curve fitting 



v(z)/vmax 



Fig. 3: Velocity profiles for the upper half of the chan- 
nel obtained by MD simulation (marker line) and fitted 
by analytical formula (solid line). In the figure, the ac- 
celeration 7^ is kept constant at 0.012 A/ps^ while the 
Knudsen number is varied. The vertical and horizontal 
axis represent normalized coordinate z/ H and velocity 
v{z)/v^a.x along the flow direction. The velocity Wmax is 
the maximal velocity found in the middle of the chan- 
nel. 



Table 2: Numerical results from MD simulations 



Kn 


7a: 
[A/ps2] 


/xx 10^ 
[Pa.s] 


[A/ps] 


Kn 


1 


CIJsHp 


fi*Kn 


A7x 


0.046 


0.0036 


1.89 


0.18 


1.87 


11.60 


21.80 




0.012 


1.91 


0.54 


1.69 


11.50 


19.40 




0.024 


2.00 


1.38 


2.20 


11.30 


25.00 


0.064 


0.0036 


1.95 


0.14 


1.97 


8.05 


15.80 




0.012 


2.02 


0.43 


1.91 


7.77 


14.90 




0.024 


2.11 


0.89 


2.10 


7.44 


15.50 


0.098 


0.0036 


2.03 


0.08 


1.90 


5.08 


9.41 




0.012 


2.31 


0.27 


2.12 


4.48 


9.51 




0.024 


2.27 


1.38 


2.26 


4.54 


10.03 



Table 3: Seidl experimental data on the couple He-Cu 
(see [S]) where TMAC depends on the collision angle of 
the gas molecules 



Angle 


Seidl 


Seidl lower 


Seidl upper 




variability 


TMAC value 


TMAC value 


10 


0.100 


0.86 


1.06 


20 


0.080 


0.81 


0.97 


30 


0.065 


0.77 


0.90 


40 


0.050 


0.73 


0.83 


50 


0.040 


0.71 


0.79 


60 


0.030 


0.69 


0.75 


70 


0.020 


0.66 


0.70 



determined. As expected, the value /3 is of order unity 
while a corresponds to an accommodation coefficient 
cr.u = 0.71 which is in the experimental range [0.66, 1.0] 
by Seidl, reported in the paper of [S] for the couple He- 
Cu. It is noted that these two bound limits 0.66 and 1.0 
correspond to the two extreme impinging at angles 70° 
and 10° of helium molecules at a copper wall in Seidl's 
experiments (see Tab. [3]). 

The impact of parameter /3 on slip velocity depends 
on the relative importance between the two terms /3 
and (Kn/i*)~^. For large Kn, the role of coefficient /3 
becomes important, e.g. up to 20% when Kn = 0.098 
while for small Kn, e.g. Kn = 0.046, this effect is almost 
negligible (6%). In order to compare with the classical 
model that involves a single parameter a, the numerical 
results are also fitted with a line passing through origin, 
the dashed line in Fig. HI The coefficient a predicted by 
the classical model takes the higher value a = 1.95 and 
shows more discrepancies with respect to the simulation 
results. 



global viscosity of the fluid /i and slip velocity VgUp are 
obtained by fitting the numerical velocity profiles (Fig. 
[3]) with equation (fT3)) . given in columns 3, 4 of Table [2J 
It is noted that the determination of /x is based the cur- 
vature of the velocity profile and the average density p, 
regardless of the redistribution of density at the steady 
state. 

From the results reported in Table [21 it can be seen 
that L*/Kn tends to increase with Kn. This trend can- 
not be accounted for by the usual one parameter model. 
In order to determine the two constants a and /3 in our 
proposed model (see Eq. [T^ , we plotted two dimension- 
less quantities '^_^'°'''' and (Kn/i*)^^ from the last two 
columns of Table |2] into Fig. SI We found the straight 
line that best fits these data. In the framework of our 
problem, the two values a = 1.81 and /3 = 0.76 were 



4 Concluding remarks and discussion 

A new slip model for flows with volume force has been 
introduced. Two parameters with physical meanings that 
relates slip velocity and velocity gradient at walls are 
suggested. The first one is the traditional tangential 
momentum accommodation coefficient, as in the widely 
used Maxwell model, while the second one accounts for 
the additional velocity that a molecule encompasses ow- 
ing to the applied force before striking at a the solid sur- 
face. Molecular dynamics calculations applied on the 
He-Cu couple allowed to determine both TMAC and 
the new parameter 13 introduced in the present work. 
TMAC is found to be in good agreement with exper- 
imental data. In the VgUp expression, the effect of pa- 
rameter j3 increases with Kn. i.e. when the limit of tran- 
sitional flow is reached. 
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Fig. 4: Relations between two dimensionless quantities 
"^-^i^lj-s. and l/(Kn/.t*). The solid line represents the equa- 
tion (fH)) which fits the numerical results. This line 
makes an angle arctana with the horizontal axis and 
cuts the vertical axis at a/3. The dashed line represents 
(fT4)) with /3 = that fits the numerical results. 



The present model is useful for the global analysis of 
flows without knowing the presence of the Knudsen 
layer. In our MD simulations, both higher density and 
slower motion of gas molecules are observed near the 
walls and cause deviation of the numerical velocity from 
the Navier-Stokes solution. Similar phenomenon was 
encountered and discussed in [5] from molecular view- 
point. Generally, when a molecule reaches a surface, it 
does not bounce back immediately but is often trapped 
by the potential well. The molecules remains near the 
wall for some time, interacts with many other solid 
atoms before escaping. Accounting for the Knudsen layer 
should leads to more accurate results (see Lockerby et 
al. |191I20| V A model involving body force as well as 
Knudsen layer may be expected in the near future. 
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